home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / rkutta.i < prev    next >
Text File  |  1995-12-12  |  12KB  |  369 lines
  1. /*
  2.    RKUTTA.I
  3.    4th order Runge-Kutta integrator (rk_integrate, rkutta)
  4.    Also Bulirsch-Stoer integrator (bs_integrate, bstoer)
  5.    After routines in Numerical Recipes by Press et.al.
  6.  
  7.    $Id$
  8.  */
  9. /*    Copyright (c) 1995.  The Regents of the University of California.
  10.                     All rights reserved.  */
  11.  
  12. /* ------------------------------------------------------------------------ */
  13.  
  14. func rk_integrate(derivative, y1, rk_x, epsilon, dx1)
  15. /* DOCUMENT y= rk_integrate(derivative, y1, x, epsilon, dx1)
  16.      integrates dydx= DERIVATIVE(y,x) beginning at (X(1),Y1) and
  17.      going to X(0) with fractional error EPSILON.  The result is
  18.      the value of y at each value in the list X.  If non-nil, DX1
  19.      will be used as initial guess for the first step size.
  20.      Otherwise, X(2)-X(1) will be the first step size guess.
  21.  
  22.      The list of X values must be monotone -- strictly increasing
  23.      or strictly decreasing; the Runge-Kutta step sizes are selected
  24.      adaptively until the next X value would be passed, when the
  25.      step size is adjusted to complete the step exactly.
  26.  
  27.      The external variable rk_maxits (default 10000) is the
  28.      maximum number of steps rk_integrate will take.
  29.  
  30.      If a function rk_yscale(y,dydx,x,dx) exists, it is used
  31.      to compute an appropriate yscale to give the EPSILON error
  32.      criterion meaning.  Otherwise, yscale is taken to be:
  33.         abs(y)+abs(dydx*dx)+1.e-30
  34.  
  35.      Based on odeint from Numerical Recipes (Press, et.al.).
  36.      If the function you are trying to integrate is very
  37.      smooth, and your X values are fairly far apart, bs_integrate
  38.      may work better than rk_integrate.
  39.  
  40.    SEE ALSO: rkutta, bs_integrate, rk_maxits, rk_minstep,
  41.              rk_maxstep, rk_ngood, rk_nbad, rkdumb, rk4
  42.  */
  43. {
  44.   if (numberof(rk_x)<2) return array(y1, dimsof(rk_x));
  45.   local rk_y;
  46.   if (is_void(dx1)) dx1= rk_x(2)-rk_x(1);
  47.   rk_nstore= -1;
  48.   rkutta, derivative, y1, rk_x(1), rk_x(0), epsilon, dx1;
  49.   return rk_y;
  50. }
  51.  
  52. func rkutta(derivative, y0,x0, x1,epsilon, dx0)
  53. /* DOCUMENT y1= rkutta(derivative, y0,x0, x1,epsilon, dx0)
  54.      integrates dydx= DERIVATIVE(y,x) beginning at (X0,Y0) and
  55.      going to X1 with fractional error EPSILON.  The result is
  56.      the value of y at X1.  DX0 will be used as the initial guess
  57.      for a step size.
  58.  
  59.      If the external variable rk_nstore is >0, rk_y and rk_x
  60.      will contain up to rk_nstore intermediate values after the
  61.      call to rkutta.  Consider using rk_integrate if you need
  62.      this feature; using rk_nstore gives you the results at
  63.      intermediate values which will tend to be closer where
  64.      the Runge-Kutta step size was smaller, while rk_integrate
  65.      forces you to specify precisely which x values you want.
  66.  
  67.      The external variable rk_maxits (default 10000) is the
  68.      maximum number of steps rkutta will take.  The variable
  69.      rk_minstep (default 0.0) is the minimum step size.  The
  70.      variable rk_maxstep (default 1.e35) is the maximum step
  71.      size, which you may need if you are storing intermediate
  72.      values (particularly with bstoer).
  73.  
  74.      If a function rk_yscale(y,dydx,x,dx) exists, it is used
  75.      to compute an appropriate yscale to give the EPSILON error
  76.      criterion meaning.  Otherwise, yscale is taken to be:
  77.         abs(y)+abs(dydx*dx)+1.e-30
  78.  
  79.      Based on odeint from Numerical Recipes (Press, et.al.).
  80.      If the function you are trying to integrate is very
  81.      smooth, bstoer will probably work better than rkutta.
  82.  
  83.    SEE ALSO: rk_integrate, bstoer, rk_nstore, rk_maxits,
  84.              rk_minstep, rk_maxstep, rk_ngood, rk_nbad, rkdumb, rk4
  85.  */
  86. {
  87.   extern rk_x, rk_y, rk_ngood, rk_nbad;
  88.   rk_ngood= rk_nbad= 0;
  89.   if (rk_nstore > 0) {
  90.     if (rk_nstore<2) rk_nstore= 2;
  91.     rk_x= array(double(x0), rk_nstore);
  92.     rk_y= array(double(y0), rk_nstore);
  93.     s= [1, 1, 1];  /* see rk_store function */
  94.   } else if (rk_nstore < 0) {
  95.     x0= rk_x(1);
  96.     x1= rk_x(0);
  97.     if (anyof(rk_x(dif)*(x1-x0) <= 0.0))
  98.       error, "given rk_x must be monotonic";
  99.     rk_y= array(double(y0), numberof(rk_x));
  100.     s= 2;
  101.   }
  102.  
  103.   dxsign= sign(x1-x0);
  104.   dx= double(abs(dx0)*dxsign);
  105.   x= double(x0);
  106.   y= double(y0);
  107.   for (n=1 ; n<=rk_maxits ; ++n) {
  108.     dydx= derivative(y, x);
  109.     if (!is_void(rk_yscale)) yscale= rk_yscale(y,dydx,x,dx);
  110.     else yscale= abs(y)+abs(dydx*dx)+1.e-30;
  111.     if (abs(dx) > rk_maxstep) dx= dxsign*rk_maxstep;
  112.     if (dxsign*(x+dx-x1) > 0.0) dx= x1-x;
  113.     if (rk_nstore<0 &&
  114.     dxsign*(x+dx-rk_x(s)) > 0.0) dx= rk_x(s)-x;
  115.     local dxdid, dxnxt;
  116.     y= rkqc(y,dydx, x,dx, derivative, epsilon,yscale, dxdid,dxnxt);
  117.     x+= dxdid;
  118.     if (dxdid == dx) ++rk_ngood;
  119.     else ++rk_nbad;
  120.     if (rk_nstore>0) s= rk_store(y,x,s);
  121.     else if (rk_nstore<0 &&
  122.          dxsign*(x-rk_x(s))>=0.0) rk_y(..,s++)= y;
  123.     all_done= (dxsign*(x-x1) >= 0.0);
  124.     if (all_done) break;
  125.     if (abs(dxnxt) < abs(rk_minstep))
  126.       error, "required step less than rk_minstep";
  127.     dx= dxnxt;
  128.   }
  129.  
  130.   if (rk_nstore>0) {
  131.     if (rk_x(s(3)) != x) {
  132.       s(2)= 1;  /* always store final value */
  133.       s= rk_store(y,x,s);
  134.     }
  135.     rk_y= rk_y(..,1:s(3));
  136.     rk_x= rk_x(1:s(3));
  137.   }
  138.   if (!all_done) error, "exceeded rk_maxits iterations";
  139.   return y;
  140. }
  141.  
  142. local rk_nstore, rk_maxits, rk_minstep, rk_maxstep, rk_ngood, rk_nbad;
  143. /* DOCUMENT rk_nstore, rk_maxits, rk_minstep, rk_maxstep,
  144.             rk_ngood, rk_nbad
  145.  
  146.      rk_nstore      maximum number of y values rkutta (bstoer) will store
  147.         after rkutta (bstoer) call, rk_y and rk_x contain stored values
  148.  
  149.      The other variables are inputs or outputs for rkutta, bstoer,
  150.      rk_integrate, or bs_integrate:
  151.  
  152.      rk_maxits      maximum number of steps (default 10000)
  153.      rk_minstep     minimum step size (default 0.0)
  154.      rk_maxstep     maximum step size (default 1.e35)
  155.      rk_ngood       number of good steps taken
  156.      rk_nbad        number of failed (but repaired) steps taken
  157.  */
  158. rk_maxits= 10000;
  159. rk_minstep= 0.0;
  160. rk_maxstep= 1.0e35;
  161. rk_nstore= 0;
  162.  
  163. func rk_store(y,x,s)
  164. {
  165.   /* s= [step number, step increment, store index] */
  166.   i= ++s(1);
  167.   if (! ((i-1)%s(2))) {
  168.     i= ++s(3);
  169.     if (i > rk_nstore) {
  170.       y2= rk_y(..,1:0:2);
  171.       x2= rk_x(1:0:2);
  172.       i= numberof(x2);
  173.       rk_y(..,1:i)= y2;
  174.       rk_x(1:i)= x2;
  175.       s(3)= ++i;
  176.       s(2)*= 2;
  177.     }
  178.     rk_y(..,i)= y;
  179.     rk_x(i)= x;
  180.   }
  181.   return s;
  182. }
  183.  
  184. func rkdumb(derivative, y0,x0, x1,nsteps, nosave=)
  185. /* DOCUMENT y_of_x= rkdumb(derivative, y0,x0, x1,nsteps)
  186.      integrates dydx= DERIVATIVE(y,x) beginning at (X0,Y0) and
  187.      going to X1 in NSTEPS 4th order Runge-Kutta steps.  The
  188.      result is dimsof(Y0)-by-(NSTEPS+1) values of y at the points
  189.      span(X0, X1, NSTEPS+1).
  190.      If the nosave= keyword is non-zero, the returned value will
  191.      simply be the final y value.
  192.  */
  193. {
  194.   dx= (x1-x0)/nsteps;
  195.   ++nsteps;
  196.   if (!nosave) y= array(double(y0), nsteps);
  197.   for (i=2 ; i<=nsteps ; ++i) {
  198.     y0= rk4(y0,derivative(y0,x0), x0,dx, derivative);
  199.     x0+= dx;
  200.     if (!nosave) y(..,i)= y0;
  201.   }
  202.   return nosave? y0 : y;
  203. }
  204.  
  205. func rkqc(y,dydx, x,dx, derivative, epsilon,yscale, &dxdid,&dxnxt)
  206. {
  207.   x0= x;
  208.   y0= y;
  209.  
  210.   for (;;) {
  211.     x= x0+dx;
  212.     if (x==x0) error, "integration step crash";
  213.     /* first take two half steps... */
  214.     dx2= 0.5*dx;
  215.     x2= x0+dx2;
  216.     y2= rk4(y0,dydx, x0,dx2, derivative);
  217.     y2= rk4(y2,derivative(y2,x2), x2,dx2, derivative);
  218.     /* ...then compare with one full step... */
  219.     y1= rk4(y0,dydx, x0,dx, derivative);
  220.     /* ...to estimate error */
  221.     y1= y2-y1;
  222.     err= max(abs(y1/yscale))/epsilon;
  223.     if (err <= 1.0) {
  224.       dxdid= dx;
  225.       dxnxt= (err>6.e-4)? 0.9*dx*err^-0.20 : 4.*dx;
  226.       break;
  227.     }
  228.     dx*= 0.9*err^-0.25;
  229.   }
  230.   return y2 + y1/15.;
  231. }
  232.  
  233. func rk4(y,dydx, x,dx, derivative)
  234. /* DOCUMENT y_at_x_plus_dx= rk4(y,dydx, x,dx, derivative)
  235.      takes a single 4th order Runge-Kutta step from X to X+DX.
  236.      DERIVATIVE(y,x) is a function returning dydx; the input DYDX
  237.      is DERIVATIVE(y,x) at the input (X,Y).  This fourth evaluation
  238.      of DERIVATIVE must be performed by the caller of rk4.
  239.  */
  240. {
  241.   dx2= 0.5*dx;
  242.   x2= x+dx2;
  243.   dydxp= derivative(y+dydx*dx2, x2);   /* slope at 1st trial midpoint */
  244.   dydxm= derivative(y+dydxp*dx2, x2);  /* slope at 2nd trial midpoint */
  245.   dydxp+= dydxm;
  246.   dydxm= derivative(y+dydxm*dx, x+dx); /* slope at trial endpoint */
  247.   return y + (dydx+dydxp+dydxp+dydxm)*(dx/6.0);
  248. }
  249.  
  250. /* ------------------------------------------------------------------------ */
  251.  
  252. func bs_integrate(derivative, y1, rk_x, epsilon, dx1)
  253. /* DOCUMENT y= bs_integrate(derivative, y1, x, epsilon, dx1)
  254.      Bulirsch-Stoer integrator, otherwise identical to rk_integrate
  255.      routine. All of the options for rk_integrate work here as well.
  256.  
  257.      Based on odeint from Numerical Recipes (Press, et.al.).
  258.      If the function you are trying to integrate is not very
  259.      smooth, or your X values are closely spaced, rk_integrate
  260.      will probably work better than bs_integrate.
  261.  
  262.    SEE ALSO: bstoer, rk_integrate, rk_maxits, rk_minstep,
  263.              rk_maxstep, rk_ngood, rk_nbad, rkdumb, rk4
  264.  */
  265. {
  266.   if (numberof(rk_x)<2) return array(y1, dimsof(rk_x));
  267.   local rk_y;
  268.   if (is_void(dx1)) dx1= rk_x(2)-rk_x(1);
  269.   rk_nstore= -1;
  270.   bstoer, derivative, y1, rk_x(1), rk_x(0), epsilon, dx1;
  271.   return rk_y;
  272. }
  273.  
  274. func bstoer(derivative, y0,x0, x1,epsilon, dx0)
  275. /* DOCUMENT y1= bstoer(derivative, y0,x0, x1,epsilon, dx0)
  276.      Bulirsch-Stoer integrator, otherwise identical to rkutta routine.
  277.      All of the options for rkutta (rk_nstore, etc.) work here as well.
  278.  
  279.      If the function you are trying to integrate is not very
  280.      smooth, rkutta will probably work better than bstoer.
  281.  
  282.    SEE ALSO: rkutta, rk_nstore, rk_maxits, rk_minstep, rk_maxstep,
  283.              rk_ngood, rk_nbad
  284.  */
  285. {
  286.   extern _rzextr_x, _rzextr_d;
  287.   rkqc= bsstep;
  288.   _rzextr_x= array(0.0, numberof(_bs_nseq));
  289.   _rzextr_d= array(double(y0), 7);
  290.   return rkutta(derivative, y0,x0, x1,epsilon, dx0);
  291. }
  292.  
  293. func bsstep(y,dydx, x,dx, derivative, epsilon,yscale, &dxdid,&dxnxt)
  294. {
  295.   x0= x;
  296.   y0= y;
  297.  
  298.   for (;;) {
  299.     for (i=1 ; i<=numberof(_bs_nseq) ; ++i) {
  300.       n= _bs_nseq(i);
  301.       y= rzextr(i, (dx/n)^2, mod_midpt(y0,dydx, x0,dx, derivative, n),
  302.         yerr, 7);
  303.       err= max(abs(yerr/yscale))/epsilon;
  304.       if (err < 1.0) {
  305.     dxdid= dx;
  306.     if (i==7) dxnxt= 0.95*dx;
  307.     else if (i==6) dxnxt= 1.2*dx;
  308.     else dxnxt= (16./*_bs_nseq(6)*/*dx)/n;
  309.     return y;
  310.       }
  311.     }
  312.     /* step failed, claimed to be unusual */
  313.     dx*= 0.0625;  /* related to numberof(_bs_nseq) */
  314.     if (x+dx == x) error, "integration step crash";
  315.   }
  316. }
  317.  
  318. _bs_nseq= [2, 4, 6, 8, 12, 16, 24, 32, 48, 64, 96];
  319.  
  320. func mod_midpt(y,dydx, x,dx, derivative, nstep)
  321. {
  322.   dx/= nstep;
  323.   ym= y;
  324.   y+= dydx*dx;
  325.   x+= dx;
  326.   dydx= derivative(y,x);
  327.   dx2= 2.*dx;
  328.   for (--nstep ; nstep ; --nstep) {
  329.     swap= ym + dydx*dx2;
  330.     ym= y;
  331.     y= swap;
  332.     x+= dx;
  333.     dydx= derivative(y,x);
  334.   }
  335.   return 0.5*(ym+y+dydx*dx);
  336. }
  337.  
  338. func rzextr(iest, xest, yest, &yerr, nuse)
  339. {
  340.   _rzextr_x(iest)= xest;
  341.   if (iest==1) {
  342.     _rzextr_d(..,1)= yest;
  343.     yerr= yest;
  344.     return yest;
  345.   } else {
  346.     m1= ((iest<nuse)? iest : nuse) - 1;
  347.     fx= _rzextr_x(iest-1 : iest-m1 : -1)/xest;  /* no more than nuse-1 */
  348.     yy= yest;
  349.     v= _rzextr_d(..,1);
  350.     _rzextr_d(..,1)= c= yy;
  351.     for (k=1 ; k<=m1 ; ++k) {
  352.       b1= fx(k)*v;
  353.       b= b1-c;
  354.       ok= double(b!=0.0);
  355.       bad= 1.0-ok;
  356.       b= ((c-v)/(b+bad))*ok;
  357.       ddy= c*b + bad*v;
  358.       c= b1*b + bad*c;
  359.       if (k!=m1) v= _rzextr_d(..,k+1);
  360.       _rzextr_d(..,k+1)= ddy;
  361.       yy+= ddy;
  362.     }
  363.     yerr= ddy;
  364.     return yy;
  365.   }
  366. }
  367.  
  368. /* ------------------------------------------------------------------------ */
  369.